Chargement du jeu de données, mise en forme et aperçu des données

Le premier jeux de données celui qui correspond au fichier CSV s’appelera “d”.

d <- read.csv("C:/Users/Cédric/Downloads/route_intersection.csv")

Transformation des dates et heures dans le bon format pour pouvoir les traiter.

d$date <- ymd_hms(d$date)

Nous faisons tout de suite une transformation de notre jeux de données pour avoir un tableau avec pour chaque date et heure toutes les informations sur les différentes intersections dans une colonne différente.

data <- d %>% mutate(
  I1 = ifelse(intersection == "I1", as.integer(vehicule), NA),
    S1 = ifelse(intersection == "S1", as.integer(vehicule), NA),
    I2 = ifelse(intersection == "I2", as.integer(vehicule), NA),
    S2 = ifelse(intersection == "S2", as.integer(vehicule), NA),
    I3 = ifelse(intersection == "I3", as.integer(vehicule), NA),
    S3 = ifelse(intersection == "S3", as.integer(vehicule), NA),
    I4 = ifelse(intersection == "I4", as.integer(vehicule), NA),
    S4 = ifelse(intersection == "S4", as.integer(vehicule), NA),
    N = ifelse(intersection == "N", as.integer(vehicule), NA)
  ) %>%
  group_by(date) %>%
  summarize(
    I1 = sum(I1, na.rm = TRUE),
    S1 = sum(S1, na.rm = TRUE),
    I2 = sum(I2, na.rm = TRUE),
    S2 = sum(S2, na.rm = TRUE),
    I3 = sum(I3, na.rm = TRUE),
    S3 = sum(S3, na.rm = TRUE),
    I4 = sum(I4, na.rm = TRUE),
    S4 = sum(S4, na.rm = TRUE),
    N = sum(N, na.rm = TRUE)
  ) %>%
  ungroup()

Calcul des données d’entrées

Calcul des données d’entrées :

t0=which(data$date=="2017-01-01 00:00:00")
n <- length(data$date)
data$E1<-data$I1+data$S2-data$I2
data$E2<-data$I2+data$S3-data$I3
data$E3 <- NA
data$E3[1:t0-1] <-data$I3[1:t0-1] - data$N[1:t0-1]
data$E3[t0:n] <- data$I3[t0:n] +data$S4[t0:n] -data$I4[t0:n] 
data$E4 <- NA
data$E4[t0:n]<-data$I4[t0:n]-data$N[t0:n]

On regarde de manière générale la circulation sur l’ensemble des intersections

data$somme <- rowSums(data[, c("E1", "E2", "E3", "E4", "N")], na.rm = TRUE)

Regroupement des données par jour

data_journaliere <- d
data_journaliere$date <- as.Date(data_journaliere$date)
data_journaliere <- data_journaliere %>% mutate(
  I1 = ifelse(intersection == "I1", as.integer(vehicule), NA),
    S1 = ifelse(intersection == "S1", as.integer(vehicule), NA),
    I2 = ifelse(intersection == "I2", as.integer(vehicule), NA),
    S2 = ifelse(intersection == "S2", as.integer(vehicule), NA),
    I3 = ifelse(intersection == "I3", as.integer(vehicule), NA),
    S3 = ifelse(intersection == "S3", as.integer(vehicule), NA),
    I4 = ifelse(intersection == "I4", as.integer(vehicule), NA),
    S4 = ifelse(intersection == "S4", as.integer(vehicule), NA),
    N = ifelse(intersection == "N", as.integer(vehicule), NA)
  ) %>%
  group_by(date) %>%
  summarize(
    I1 = max(I1, na.rm = TRUE),
    S1 = max(S1, na.rm = TRUE),
    I2 = max(I2, na.rm = TRUE),
    S2 = max(S2, na.rm = TRUE),
    I3 = max(I3, na.rm = TRUE),
    S3 = max(S3, na.rm = TRUE),
    I4 = ifelse(all(is.na(I4)), 0, max(I4, na.rm = TRUE)),
    S4 = ifelse(all(is.na(S4)), 0, max(S4, na.rm = TRUE)),
    N = max(N, na.rm = TRUE)
  ) %>%
  ungroup()

t0_jour=which(data_journaliere$date=="2017-01-01")
n <- length(data_journaliere$date)
data_journaliere$E1 <- data_journaliere$I1 + data_journaliere$S2 - data_journaliere$I2
data_journaliere$E2 <- data_journaliere$I2 + data_journaliere$S3 - data_journaliere$I3
data_journaliere$E3 <- NA
data_journaliere$E3[1:t0_jour-1] <- data_journaliere$I3[1:t0_jour-1] - data_journaliere$N[1:t0_jour-1]
data_journaliere$E3[t0_jour:n] <- data_journaliere$I3[t0_jour:n] + data_journaliere$S4[t0_jour:n] - data_journaliere$I4[t0_jour:n] 
data_journaliere$E4 <- NA
data_journaliere$E4[t0_jour:n] <- data_journaliere$I4[t0_jour:n]-data_journaliere$N[t0_jour:n]
data_journaliere$somme <- rowSums(data_journaliere[, c("E1", "E2", "E3", "E4", "N")], na.rm = TRUE)

Regroupement des données par mois

data_mois = data_journaliere %>% 
  mutate(month = format(date, "%Y-%m")) %>%
  group_by(month) %>%
  summarize(
    I1 = sum(I1, na.rm = TRUE),
    S1 = sum(S1, na.rm = TRUE),
    I2 = sum(I2, na.rm = TRUE),
    S2 = sum(S2, na.rm = TRUE),
    I3 = sum(I3, na.rm = TRUE),
    S3 = sum(S3, na.rm = TRUE),
    I4 = sum(I4, na.rm = TRUE),
    S4 = sum(S4, na.rm = TRUE),
    N = sum(N, na.rm = TRUE),
    E1 = sum(E1),
    E2 = sum(E2),
    E3 = sum(E3),
    E4 = sum(E4,na.rm=T), 
    somme = sum(somme, na.rm = TRUE)
  )

Observation générale des données

ggplot(data, aes(date,somme)) +
  geom_point()

On peut déjà remarquer une augmentation de la circulation de manière générale sur l’ensemble des intersections.

summary(data_journaliere)
##       date                  I1               S1              I2      
##  Min.   :2015-11-01   Min.   : 18.00   Min.   : 3.00   Min.   : 5.0  
##  1st Qu.:2016-03-31   1st Qu.: 42.00   1st Qu.: 8.00   1st Qu.:15.0  
##  Median :2016-08-30   Median : 62.00   Median :12.00   Median :19.0  
##  Mean   :2016-08-30   Mean   : 63.35   Mean   :12.43   Mean   :20.9  
##  3rd Qu.:2017-01-29   3rd Qu.: 81.00   3rd Qu.:16.00   3rd Qu.:24.0  
##  Max.   :2017-06-30   Max.   :156.00   Max.   :34.00   Max.   :48.0  
##                                                                      
##        S2               I3              S3              I4        
##  Min.   : 2.000   Min.   :  6.0   Min.   :  3.0   Min.   : 0.000  
##  1st Qu.: 4.000   1st Qu.: 16.0   1st Qu.:  8.0   1st Qu.: 0.000  
##  Median : 6.000   Median : 22.0   Median : 12.0   Median : 0.000  
##  Mean   : 6.441   Mean   : 28.5   Mean   : 19.3   Mean   : 5.362  
##  3rd Qu.: 7.250   3rd Qu.: 33.0   3rd Qu.: 21.0   3rd Qu.:15.000  
##  Max.   :19.000   Max.   :180.0   Max.   :164.0   Max.   :36.000  
##                                                                   
##        S4               N                E1               E2       
##  Min.   : 0.000   Min.   : 0.000   Min.   : 11.00   Min.   : 3.00  
##  1st Qu.: 0.000   1st Qu.: 6.000   1st Qu.: 32.00   1st Qu.: 8.00  
##  Median : 0.000   Median : 8.000   Median : 48.00   Median :10.00  
##  Mean   : 2.053   Mean   : 9.523   Mean   : 48.89   Mean   :11.69  
##  3rd Qu.: 4.000   3rd Qu.:10.000   3rd Qu.: 64.00   3rd Qu.:14.00  
##  Max.   :18.000   Max.   :62.000   Max.   :138.00   Max.   :36.00  
##                                                                    
##        E3               E4            somme       
##  Min.   :  3.00   Min.   : 0.00   Min.   : 26.00  
##  1st Qu.:  9.00   1st Qu.: 9.00   1st Qu.: 62.00  
##  Median : 13.00   Median :10.00   Median : 86.00  
##  Mean   : 17.99   Mean   :10.23   Mean   : 91.14  
##  3rd Qu.: 20.00   3rd Qu.:12.00   3rd Qu.:117.00  
##  Max.   :172.00   Max.   :24.00   Max.   :318.00  
##                   NA's   :427
mean(data_journaliere$I1)>sum(mean(data_journaliere$I2)+mean(data_journaliere$I3)+mean(data_journaliere$I4))
## [1] TRUE

Si l’on observe bien les colonnes on peut voir que I1 est l’intersection avec le nombre moyen de véhicule par jour le plus élevé, ce nombre est même supérieur à la somme des moyennes (journalières) de I2,I3,I4. Nous en reparlerons dans le rapport …

Ici nous avons réalisé un modèle multiplicatif à l’aide de la fonction decompose.

On a un ggplot avec des résidus en couleur et un modèle qui suit les données assez bien.

PS: les résidus semblent incorrects

## [1] NA
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

Ici nous avons cherché à observer la tendance globale sur les entrées (E1, E2, E3) avec un simple modèle linéaire :

## `geom_smooth()` using formula = 'y ~ x'

Jours de la semaine / Week-end.

Nous avons eu l’intuition de regarder différemment les jours de la semaine et les week-ends parce qu’en observant plus précisemment les données on se rend compte qu’il y a moins de circulation pendant les week-end, alors pour mettre en place une période il nous a semblé pertinent de séparer les deux.

Voici le graphique qui nous a permis de déduire cela :

plot_ly(data[0:500,], x = ~date, y = ~somme, type = 'scatter', mode = 'markers')

On temarque bien que les données pour les week-end sont plus basses.

data$jour <- weekdays(data$date)
week_end <- which(data$jour == "samedi" | data$jour == "dimanche")  
data$week_end <- FALSE
data$week_end[week_end] <- TRUE

Aperçu avec un graphique.

ggplot(data, aes(date, somme, color = week_end)) +
  geom_line() +
  ggtitle("Données en fonction des jours de la semaine")
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

On observe bien la différence, on voit bien que les valeurs sont plus faibles pour les week-end. Cependant en reprenant le dernier graphique nous voyons bien qu’il y a une période beaucoup plus marquée par semaine que par jour.

Données de l’intersection 1

ggplot(data = data, aes(date, I1)) +
  geom_line(aes(color = "I1"), linetype = "solid") +
  geom_line(aes(y = E1, color = "E1"), linetype = "solid") +
  geom_line(aes(y = S1, color = "S1"), linetype = "solid") +
  scale_color_manual(values = c("I1" = "lightskyblue3", "E1" = "indianred", "S1" = "khaki3" )) +
  ggtitle("Intersection 1")

On remarque que les données de I1 sont plus élevé que celle de E1 mais que cela suit la même tendance. On peut aussi remarquer qu’il y a beaucoup plus d’entrées que de sorties dans cette intersection.

Donnnées de l’intersection 2

ggplot(data = data, aes(date,I2)) +
  geom_line(aes(color = "I2")) +
  geom_line(aes(y=E2, color ="E2")) +
  geom_line(aes(y = S2, color = "S2")) +
  scale_color_manual(values = c("I2" = "lightskyblue3", "E2" = "indianred", "S2" = "khaki3" )) +
  ggtitle("Intersection 2")

Sur ce graphique on observe beaucoup mieux la cassure et de même que le graphique pour l’intersection 1 on remarque que peut importe si c’est l’entrée, la sortie, les données suivent le même modèle.

Données de l’intersection 3

ggplot(data = data, aes(date,I3)) +
  geom_line(aes(color = "I3")) +
  geom_line(aes(y=E3, color ="E3")) +
  geom_line(aes(y = S3, color = "S3")) +
  scale_color_manual(values = c("I3" = "lightskyblue3", "E3" = "indianred", "S3" = "khaki3" )) +
  ggtitle("Intersection 3")

Il n’ya a pas vraiment de tendance globale pour ce graphique, c’est une osillation perpétuelle peut-être peut-on trouver une périodicié mais cela ne semble par régulier. On pourra quand même calculer la tendance linéaire pour avoir une tendance générale, plutôt d’augmentation ou de diminution.

Données de l’intersection 4 (à partir du 1er janvier 2017)

ggplot(data = data[t0:length(data$date),], aes(date,I4)) +
  geom_line(aes(color = "I4")) +
  geom_line(aes(y=E4, color ="E4")) +
  geom_line(aes(y = S4, color = "S4")) +
  scale_color_manual(values = c("I4" = "lightskyblue3", "E4" = "indianred", "S4" = "khaki3" )) +
  ggtitle("Intersection 4")

On observe une augmentation progressive. Puis cela semble stagner.

Lissage et modélisation sur les données par jour de chaque intersection :

Intersection 1

Tests série différence, normalisation pas très intéressant

Moyenne mobile

Nous avons choisi un k grand parce qu’il y a beaucoup de données. La barre verticale bleue correspond à l’ouverture de la 4e intersection.

k <- 50
res_rollmean_E1 <- rollmean(data$E1,k = 2*k+1)
res_rollmean_E1 <- c(rep(NA,k),res_rollmean_E1,rep(NA,k))
ggplot(data,aes(x=date,y=E1)) +
  geom_point(size = 0.5, col = "grey") +
  geom_line(aes(y=res_rollmean_E1),col="red",size=0.5) + 
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Moyenne mobile pour l'intersection 1")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.
## Warning: Removed 100 rows containing missing values (`geom_line()`).

On peut observer une tendance à la hausse jusqu’au premier janvier 2017 et ensuite les données semblent plutôt stagner.

Lissage exponentielle simple :

res_ses_E1 <- ses(data$E1,alpha = .1,initial = "simple")
res_lissage_exp_E1 <- res_ses_E1$fitted
ggplot(data, aes(date, E1)) +
  geom_point(size = 0.5, col = "grey") +
  geom_line(aes(y = res_lissage_exp_E1), size = 0.5, col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Lissage exponentiel simple pour l'intersection 1")
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.

Le lissage exponentiel ne nous apprend rien de plus. On peut essayer d’en faire un double.

Lissage exponentiel double :

res_ses_double_E1 <- ses(res_lissage_exp_E1,alpha = .1,initial = "simple")
res_lissage_exp_double_E1 <- res_ses_double_E1$fitted
ggplot(data, aes(date, E1)) +
  geom_point(size = 0.5, col = "grey") +
  geom_line(aes(y = res_lissage_exp_double_E1), size = 0.5, col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Lissage exponentiel double pour l'intersection 1")
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.

On voit peut-être mieux le côté étalement des données.

Pour les modélisations nous avons essayés le modèle multiplicatif parce que c’est celui qui ressemble le mieux aux données : Nous avons choisi une fréquence de 24 pour les 24 heures d’une journée.

Tendance linéaire

Essayons de modéliser une tendance linéaire pour avoir un aperçu génréral de l’augmentation des données.

attach(data)
## L'objet suivant est masqué _par_ .GlobalEnv:
## 
##     week_end
res_lm_E1 <- lm(E1~date)
ggplot(data, aes(date,E1)) +
  geom_point(size = 0.5,col = "grey") +
  geom_line(aes(y = res_lm_E1$fitted.values), col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Tendance linéaire sur l'intersection 1")
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.

On le voyait déjà avant mais les chiffres parlent bien il y a bien une augmentation générale des données.

Tendance linéaire avant et après l’ouverture de la 4e intersection :

Pour cela séparons le jeux de données :

n <- nrow(data)
data_avant <- data[0:t0,]
data_apres <- data[(t0+1):n,]

mettons en place les deux modèles

Valeur dans ‘data’ :

data$tend_E1 <- NA
data$tend_E1[0:t0] <- data_avant$tend_E1
data$tend_E1[(t0+1):n] <- data_apres$tend_E1

Visualisation

ggplot(data, aes(date, E1)) +
  geom_point(col = "grey") +
  geom_line(aes(y = data$tend_E1), col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Tendance linéaire avant et après l'ouverture de la 4e intersection")
## Warning: Use of `data$tend_E1` is discouraged.
## ℹ Use `tend_E1` instead.
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.

On observe bien une différence d’augmentation des données.

Modèle multiplicatif.

data$E1_ts <- ts(data$E1,frequency = 24)
res_decompose_E1 <- decompose(data$E1_ts,type="multiplicative")
data$fitted <- as.numeric(res_decompose_E1$trend*res_decompose_E1$seasonal)
ggplot(data = data, aes(date, E1))+
  geom_line()+
  geom_line(aes(y = fitted), col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Modèle multiplicatif pour l'intersection 1")
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.
## Warning: Removed 24 rows containing missing values (`geom_line()`).

On voit bien ici que la modélisation avec un modèle multiplicatif correspond bien aux données mais la fonction ‘decompose’ de R se base sur la moyenne mobile donc on n’obtient pas un modèle très propre mais plutôt un lissage. Cependant ce modèle s’adapte très bien aux données.

Nous avons aussi essayé de créer nous même un modèle multiplicatif :

modele_multi <- function(serie_vehicule,serie_date,p){
  logar=lm(log(serie_vehicule)~serie_date) 
  detendance=logar$residuals
  logar_fitted=logar$fitted.values
  
  n <- length(detendance)
  
  N <- floor(n/p)
  
  periode <- rep(NA,p)
  
  for(i in 1:p){
    # Calculer la moyenne de chacune des étapes de l’ensemble des périodes
    periode[i] <- mean(detendance[i+(1:(N+1)-1)*p],na.rm=T)
  }
  if(n == N*p) periode <- rep(periode,N)
  
  if(n > N*p) periode <- c(rep(periode,N),periode[1:(n-(N*p))])

  modele <- exp(periode) * exp(logar_fitted)
  
  residu <- serie_vehicule - modele #problème avec résidu
  
  tableau <- data.frame(modele, residu)
  return(tableau)
}

Nous l’avons tester sur les données :

mod_multiplicatif_E1 <- modele_multi(data$E1, data$date,24)
ggplot(data, aes(date, E1)) +
  geom_point(size = 0.5, col = "grey") +
  geom_line(aes(y = mod_multiplicatif_E1$modele), col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Autre modèle multiplicatif")
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.

Le modèle correspond plutôt bien aux données, mais nous n’avons un problème avec les résidus. C’es pour cela que dans la suite nous avons préféré nous concentrer sur les données par jour plutôt que par heure.

Intersection 2 :

ggplot(data, aes(date, E2)) + geom_point(size = 0.5, col = "grey") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Données de l'intersection 2")
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.

Les données parraissent étonantes dans un premier temps mais c’est en fait parce que l’échelle est petite, il y peut de voitures en moyenne qui passent par l’intersection 2.

Moyenne mobile

k <- 50
res_rollmean_E2 <- rollmean(data$E2,k = 2*k+1)
res_rollmean_E2 <- c(rep(NA,k),res_rollmean_E2,rep(NA,k))
ggplot(data,aes(x=date,y=E2)) +
  geom_point(size = 0.5, col = "grey") +
  geom_line(aes(y=res_rollmean_E2),col="red",size=0.5) + 
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Moyenne mobile sur l'intersection 1")
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.
## Warning: Removed 100 rows containing missing values (`geom_line()`).

On peut ainsi observer la tendance globale, on voit très nettement une cassure dans les données ici.

Lissage exponentiel simple

res_ses_E2 <- ses(data$E2,alpha = .1,initial = "simple")
res_lissage_exp_E2 <- res_ses_E2$fitted
ggplot(data, aes(date, E2)) +
  geom_point(size = 0.5, col = "grey") +
  geom_line(aes(y = res_lissage_exp_E2), size = 0.5, col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Lissage exponentiel sur les données E2")
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.

Cela nous apprend pas grand chose de plus. La moyenne mobile semble un peu plus significative. On peut jsute dire qu’il y a un étalement des données après la date de l’ouverture de la dernière intersection.

Tendance linéaire

Comme pour l’intersection I1 nous allons regarder la tendance linéaire.

attach(data)
## L'objet suivant est masqué _par_ .GlobalEnv:
## 
##     week_end
## Les objets suivants sont masqués depuis data_apres:
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, week_end
## Les objets suivants sont masqués depuis data_avant:
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, week_end
## Les objets suivants sont masqués depuis data (pos = 5):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, week_end
res_lm_E2 <- lm(E2~date)
ggplot(data, aes(date,E2)) +
  geom_point(col = "grey") +
  geom_line(aes(y = res_lm_E2$fitted.values), col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Tendance linéaire intersection 2")
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.

Contrairement à l’intersection 1 on se rend bien compte que ce modèle modélise très mal les données.

plot(res_lm_E2$residuals, type ="p")
abline(h=0, col = "red")

On remarque bien qu’il y a un moment où les résidus passent tous en dessous de 0, ce qui est le signe d’une cassure.

Essayons de modéliser séparement avant et après 1er janvier 2017.

Valeur dans ‘data’ :

data$tend_E2 <- NA
data$tend_E2[0:t0] <- data_avant$tend_E2
data$tend_E2[(t0+1):n] <- data_apres$tend_E2

Visualisation

ggplot(data, aes(date, E2)) +
  geom_point(col = "grey") +
  geom_line(aes(y = data$tend_E2), col = "red")
## Warning: Use of `data$tend_E2` is discouraged.
## ℹ Use `tend_E2` instead.

On voit très bien la cassure ici. Nous avons donc une augmentation plus importante de la circulation après la mise en place de la nouvelle intersection qu’avant. Mais avec un étalement des données, une plus grande répartition.

Intersection 3

En observant les données de l’intersection 3 on voit déjà une échelle beaucoup plus grande et un graphique qui ne ressemble pas du tout aux précédents.

Moyenne mobile

k <- 50
res_rollmean_E3 <- rollmean(data$E3,k = 2*k+1)
res_rollmean_E3 <- c(rep(NA,k),res_rollmean_E3,rep(NA,k))
ggplot(data,aes(x=date,y=E3)) +
  geom_point(size = 0.5, col = "grey") +
  geom_line(aes(y=res_rollmean_E3),col="red",size=0.5) + 
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Moyenne mobile pour l'intersection 3")
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.
## Warning: Removed 100 rows containing missing values (`geom_line()`).

Tendance qui a plutôt l’air de stagner mais avec un très fort déséquilibre entre les données. De plus nous n’avons pas l’impression d’observer un modèle qui se répète il parait difficile de mettre en place un périodicité.

Lissage exponentiel simple

Peut-être que le lissage exponentiel peut nous en apprendre un peu plus.

res_ses_E3 <- ses(data$E3,alpha = .1,initial = "simple")
res_lissage_exp_E3 <- res_ses_E3$fitted
ggplot(data, aes(date, E3)) +
  geom_point(size = 0.5, col = "grey") +
  geom_line(aes(y = res_lissage_exp_E3), size = 0.5, col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Lissage exponentiel pour l'intersection 3")
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.

Il n’y a pas beaucoup de différence.

Tendance linéaire

Observons la tendance linéaire pour visualiser de manière générale une potentielle augmentation ou diminution de la circulation :

attach(data)
## L'objet suivant est masqué _par_ .GlobalEnv:
## 
##     week_end
## Les objets suivants sont masqués depuis data_apres (pos = 3):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, tend_E1, week_end
## Les objets suivants sont masqués depuis data_avant (pos = 4):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, tend_E1, week_end
## Les objets suivants sont masqués depuis data (pos = 5):
## 
##     anomalie, date, E1, E1_ts, E2, E3, E4, fitted, I1, I2, I3, I4,
##     jour, N, residu, S1, S2, S3, S4, somme, tend_E1, week_end
## Les objets suivants sont masqués depuis data_apres (pos = 6):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, week_end
## Les objets suivants sont masqués depuis data_avant (pos = 7):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, week_end
## Les objets suivants sont masqués depuis data (pos = 8):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, week_end
res_lm_E3 <- lm(E3~date)
ggplot(data, aes(date,E3)) +
  geom_point(col = "grey", size = 0.5) +
  geom_line(aes(y = res_lm_E3$fitted.values), col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Tendance linéaire pour l'intersection 3")
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.

Il semble y avoir une légère augmentation générale. Comme précédemment on peut essayer de regarder avant et après la mise en place de la dernière intersection.

ggplot(data, aes(date, E3)) +
  geom_point(col = "grey", size = 0.5) +
  geom_line(aes(y = data$tend_E3), col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Tendance linéaire avant et après l'ouverture de l'intersection 4 pour les données de la 3e")
## Warning: Use of `data$tend_E3` is discouraged.
## ℹ Use `tend_E3` instead.
## Warning: Use of `data$date` is discouraged.
## ℹ Use `date` instead.

Alors c’est très intéressant parce qu’après la mise en place de la 4e intersection les données semblent plutôt diminuer !

Analyse des résidus :

Analysons les résidus de ces deux tendances tendances linéaires pour voir si l’on observe pas quelques chose :

plot(res_lm_1_E3$residuals)

plot(res_lm_2_E3$residuals)

Rien de réellement intéressant, nous verrons avec les données par jour.

Intersection I4

Moyenne mobile

k <- 50
res_rollmean_E4 <- rollmean(data$E4[t0:n],k = 2*k+1)
res_rollmean_E4 <- c(rep(NA,k),res_rollmean_E4,rep(NA,k))
ggplot(data[t0:n,],aes(x=date,y=E4)) +
  geom_point(size = 0.5, col = "grey") +
  geom_line(aes(y=res_rollmean_E4),col="red",size=0.5) + 
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Moyenne mobile pour l'intersection 4")
## Warning: Removed 100 rows containing missing values (`geom_line()`).

On observe bien une légère augmentation au début et puis les données commences à stager.

Lissage exponentiel simple

res_ses_E4 <- ses(data$E4[t0:n],alpha = .1,initial = "simple")
res_lissage_exp_E4 <- res_ses_E4$fitted
ggplot(data[t0:n,], aes(date, E4)) +
  geom_point(size = 0.5, col = "grey") +
  geom_line(aes(y = res_lissage_exp_E4), size = 0.5, col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("lissage exponentiel pour l'intersection 4")

Peu d’informations en plus par rapport à la moyenne mobile.

Tendance linéaire

attach(data[t0:n,])
## L'objet suivant est masqué _par_ .GlobalEnv:
## 
##     week_end
## Les objets suivants sont masqués depuis data_apres (pos = 3):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, tend_E1, tend_E2, week_end
## Les objets suivants sont masqués depuis data_avant (pos = 4):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, tend_E1, tend_E2, week_end
## Les objets suivants sont masqués depuis data (pos = 5):
## 
##     anomalie, date, E1, E1_ts, E2, E3, E4, fitted, I1, I2, I3, I4,
##     jour, N, residu, S1, S2, S3, S4, somme, tend_E1, tend_E2, week_end
## Les objets suivants sont masqués depuis data_apres (pos = 6):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, tend_E1, week_end
## Les objets suivants sont masqués depuis data_avant (pos = 7):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, tend_E1, week_end
## Les objets suivants sont masqués depuis data (pos = 8):
## 
##     anomalie, date, E1, E1_ts, E2, E3, E4, fitted, I1, I2, I3, I4,
##     jour, N, residu, S1, S2, S3, S4, somme, tend_E1, week_end
## Les objets suivants sont masqués depuis data_apres (pos = 9):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, week_end
## Les objets suivants sont masqués depuis data_avant (pos = 10):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, week_end
## Les objets suivants sont masqués depuis data (pos = 11):
## 
##     anomalie, date, E1, E2, E3, E4, fitted, I1, I2, I3, I4, jour, N,
##     residu, S1, S2, S3, S4, somme, week_end
res_lm_E4 <- lm(E4~date)
ggplot(data[t0:n,], aes(date,E4)) +
  geom_point(col = "grey", size = 0.5) +
  geom_line(aes(y = res_lm_E4$fitted.values), col = "red") +
  geom_line(aes(x=data$date[t0]), col = "blue") +
  ggtitle("Tendance linéaire pour l'intersection 4")

Conclusion

Nous avons eu un aperçu général des données, nous avons commencé à voir qu’il y a certainement une cassure au niveau de l’apparition de la dernière intersection pour les trois premières intersections. Pour continuer cette analyse afin de la rendre plus facile nous avons travaillé sur les données par jour afin de pouvoir dégager une potentiel période de 7 jours et pour avoir un meilleur aperçu avec moins de données.

Graphique sur l’ensemble des données

La moyenne mobile par jour

Création de colonnes pour mettre la moyenne mobile

k <- 3
data_journaliere$rmean_E1 <- rollmean(data_journaliere$E1, k = 2*k+1, fill = NA)
data_journaliere$rmean_E2 <- rollmean(data_journaliere$E2, k = 2*k+1, fill = NA)
data_journaliere$rmean_E3 <- rollmean(data_journaliere$E3, k = 2*k+1, fill = NA)
data_journaliere$rmean_somme =  rollmean(data_journaliere$somme, k = 2*k+1, fill = NA)

Visualisation : Avec ou sans points :

ggplot(data=data_journaliere, aes(x=date)) +
  geom_line(aes(y=rmean_E1, color="E1"), size=1) +
  geom_line(aes(y=rmean_E2, color="E2"), size=1) +
  geom_line(aes(y=rmean_E3, color="E3"), size=1) +
  geom_point(aes(y=E1, color="E1"), size=0.2, alpha=0.5) +
  geom_point(aes(y=E2, color="E2"), size=0.2, alpha=0.5) +
  geom_point(aes(y=E3, color="E3"), size=0.2, alpha=0.5) +
  labs(title="Moyennes Mobiles (k=10) et Points des Variables E1, E2 et E3",
       x="Date",
       y="Valeur") +
  scale_color_manual(values=c("E1"="blue", "E2"="green", "E3"="red")) +
  theme_minimal()+
  ylim(0,100)+
  geom_vline(xintercept=as.Date("2017-01-01"), linetype="dashed", color="black")
## Warning: Removed 6 rows containing missing values (`geom_line()`).
## Removed 6 rows containing missing values (`geom_line()`).
## Removed 6 rows containing missing values (`geom_line()`).
## Warning: Removed 8 rows containing missing values (`geom_point()`).
## Warning: Removed 3 rows containing missing values (`geom_point()`).

ggplot(data=data_journaliere, aes(x=date)) +
  geom_line(aes(y=rmean_E1, color="E1"), size=1) +
  geom_line(aes(y=rmean_E2, color="E2"), size=1) +
  geom_line(aes(y=rmean_E3, color="E3"), size=1) +
  labs(title="Moyennes Mobiles (k=10) et Points des Variables E1, E2 et E3",
       x="Date",
       y="Valeur") +
  scale_color_manual(values=c("E1"="blue", "E2"="green", "E3"="red")) +
  theme_minimal()+
  ylim(0,100)+
  geom_vline(xintercept=as.Date("2017-01-01"), linetype="dashed", color="black")
## Warning: Removed 6 rows containing missing values (`geom_line()`).
## Warning: Removed 6 rows containing missing values (`geom_line()`).
## Removed 6 rows containing missing values (`geom_line()`).

On remarque que c’est l’entrée 1 qui est la plus engorgée.

Sur l’ensemble des données la colonne ‘somme’ :

ggplot(data=data_journaliere,aes(x=date))+
  geom_point(aes(y=somme),color="purple",size=0.9)+
  geom_line(aes(y=rmean_somme),size=0.6)+
  ylim(0,150)+ labs(title="Moyennes Mobiles (k=3) sur somme",
       x="Date",
       y="Valeur")
## Warning: Removed 39 rows containing missing values (`geom_point()`).
## Warning: Removed 6 rows containing missing values (`geom_line()`).

Modèle multiplicatif

Ici on exploite la fonction modele_multi afin de calculer un modèle multiplicatif puis on le représente graphiquement.

La période tester ici est de 24 heures sur les données par heure donc data

Ici la meme chose mais le graphique à une échelle “dézoomer” sur toute la période étudier

## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

Nous avons tenté un modèle linéaire sur l’ensemble des données, c’est une régression linéaire sur les résidus d’une regression linéaire. L’idée ici était de voir si en enlevant la tendance globale d’augmentation de circulation on pouvait observer une baisse (relative) liée à l’introduction de la nouvelle intersection I4. Finalement l’approche ne nous a pas semblé très concluante…

## `geom_smooth()` using formula = 'y ~ x'

Modélisation pour I1

Modèle additif

En regroupant les données par jour cela nous permet de modéliser les données E1 avec un modèle additif, étant donné que le modèle multiplicatif ne nous permet pas de faire des prévision.

calcul_periodicite <- function(serie,p){
n <- length(serie)
N <- floor(n/p)
periode <- rep(NA,p)
for(i in 1:p){
periode[i] <- mean(serie[i+(1:(N+1)-1)*p],na.rm=T)
}
if(n == N*p) periode <- rep(periode,N)
if(n > N*p) periode <- c(rep(periode,N),periode[1:(n-(N*p))])
return(periode)
}
res_lm_E1 <- lm(data_journaliere$E1~data_journaliere$date)
E1_detendance <- res_lm_E1$residuals 
periode_E1 <- calcul_periodicite(E1_detendance, 7)
data_journaliere$modele_additif_E1 <- res_lm_E1$fitted.values + periode_E1

Visualisation :

ggplot(data_journaliere, aes(date, E1)) +
  geom_point(size = 0.5) +
  geom_line(aes(y = modele_additif_E1), col = "red")

On peut voir que ça modélise plutôt pas trop mal les données. Analysons maintenant les résidus pour voir les anomalie.

Anomalies

residus <- data_journaliere$E1 - data_journaliere$modele_additif_E1

Regardons maintenant par rapport aux quantiles :

q1 = quantile(probs = 0.99, residus, na.rm = TRUE)
index=which(residus > q1)
data_journaliere$anomalie <- ifelse(data_journaliere$E1 %in% data_journaliere$E1[index], TRUE, FALSE)

 ggplot(data_journaliere,aes(x = date,y = E1,color=anomalie))+
   geom_point()+
   geom_line(aes(y=modele_additif_E1), col = "red") +
   geom_vline(xintercept=as.Date("2017-01-01"), linetype="dashed", color="black")+
  labs(title="Modèle additif sur une période de 7 jours",
       x="Date",
       y="Nombre véhicules")

On observe qu’il y a plus d’anomalies après l’ouverture d’I4.

Cassure test de Chow

Nous avons cherché à confirmer l’hypothèse que les données modéliser précedemment présente une cassure “significative”. On utilisera la fonction chow.test() du package gap.

Pour cela nous avons réalisé un test de chow en séparant les données avant l’introduction d’I4 : soit avant t0_jour = “2017-01-01” et les données après cette date. En utilisant encore un modèle additif (période= 7j)…

Le test est concluant puisque la statistique de test est non comprise dans la région de test adpatée à nos données (degré de liberté 1 et 2) avec un risque de 5% (de quantile fisher(paramètres) à +Inf)

t0_jour=which(data_journaliere$date=="2017-01-01")
n <- length(data_journaliere$date)

tendance_E1_1=lm(data_journaliere$E1[1:t0_jour]~data_journaliere$date[1:t0_jour])$fitted
tendance_E1_2=lm(data_journaliere$E1[t0_jour:n]~data_journaliere$date[t0_jour:n])$fitted


res_chow<- chow.test(tendance_E1_1,data_journaliere$date[1:t0_jour],tendance_E1_2,data_journaliere$date[t0_jour:n])
k=2
f=c(qf(0.95,df1=k,df2=(length(tendance_E1_2)+length(tendance_E1_1)+2*k)),Inf)
if (res_chow[1] > f[1]){
  print(TRUE)
} else {
  print(FALSE)
}
## [1] FALSE

La cassure n’étant confirmé par le test ci-dessus, nous avons fait des prédictions.

Prédictions

periode_E1[1:7]
## [1] -12.984780   8.792783   5.305978   5.485840   4.918575   1.512230 -13.182145
periode_E1[length(periode_E1)]
## [1] 1.51223
#periode j6
dates_futures <- data.frame(date = seq(max(data_journaliere$date) + 1,
                                       by = 1, length.out = 40))

Pour prédire des valeurs futures, vous devez avoir les dates futures disponibles. Créons un nouveau dataframe avec les nouvelles dates pour la prédiction. Supposons que vous voulez prédire les 40 prochains jours.

dates_futures <- data.frame(date = seq(max(data_journaliere$date) + 1,
                                       by = 1, length.out = 40))

Prédictions :

modele <- lm(E1 ~ date, data = data_journaliere)
predictions <- predict(modele, newdata = dates_futures)
predictions_E1=rep(NA,40)
for (i in 1:40){
  valeur_periode=5+i%%7
  predictions_E1[i] = predictions[i] + periode_E1[valeur_periode]

}
prediction_E1 = predictions + c(periode_E1[6:7],periode_E1[1:7],periode_E1[1]) #sur 10

derniere_date <- max(data_journaliere$date, na.rm = TRUE)

Créer une séquence de 10 nouvelles dates à partir du jour après la dernière date

nouvelles_dates <- seq(derniere_date + 1, by = "days", length.out = 40)

Créer un nouveau dataframe avec les nouvelles dates

nouvelles_lignes <- data.frame(date = nouvelles_dates)

Ajouter des NA aux autres colonnes du nouveau dataframe

for(col in names(data_journaliere)) {
  if(col != "date") {
    nouvelles_lignes[[col]] <- NA
  }
}

Combiner le dataframe original avec le nouveau dataframe de nouvelles dates

data_journaliere2 <- rbind(data_journaliere, nouvelles_lignes)
ggplot(data_journaliere2, aes(x = date, y = E1)) +
  geom_smooth(method="lm")+
  geom_point(aes(color = anomalie)) +
  scale_color_manual(values = c("TRUE" = "red", "FALSE" = "blue")) +
  geom_line(aes(y = c(data_journaliere$modele_additif_E1, predictions_E1)), col = "grey40") +
  geom_vline(xintercept = as.Date("2017-01-01"), linetype = "dashed", color = "black") +
  labs(title = "Modèle additif sur une période de 7 jours (E1)",
       x = "Date",
       y = "Nombre véhicules") 
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 40 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 40 rows containing missing values (`geom_point()`).

Intersection 2

Modèle additif

res_lm_E2 <- lm(data_journaliere$E2~data_journaliere$date)
E2_detendance <- res_lm_E2$residuals 
periode_E2 <- calcul_periodicite(E2_detendance, 7)
data_journaliere$modele_additif_E2 <- res_lm_E2$fitted.values + periode_E2
ggplot(data_journaliere, aes(date, E2)) +
  geom_point(size = 0.5) +
  geom_line(aes(y = modele_additif_E2), col = "red")

Mauvais modèle :

cassure modèle linéaire :

t0_jour=which(data_journaliere$date=="2017-01-01")
n <- length(data_journaliere$date)

tendance_E2_1=lm(data_journaliere$E2[1:t0_jour]~data_journaliere$date[1:t0_jour])$fitted
tendance_E2_2=lm(data_journaliere$E2[t0_jour:n]~data_journaliere$date[t0_jour:n])$fitted


res_chow_2<- chow.test(tendance_E2_1,data_journaliere$date[1:t0_jour],tendance_E2_2,data_journaliere$date[t0_jour:n])
k=2
f=c(qf(0.95,df1=k,df2=(length(tendance_E2_2)+length(tendance_E1_1)+2*k)),Inf)
if (res_chow_2[1] > f[1]){
  print(TRUE)
} else {
  print(FALSE)
}
## [1] TRUE

Il y a donc une cassure bien visible sur l’intersection 2.

visualisation :

data_journaliere$tendance_lin_E2 <- NA
data_journaliere$tendance_lin_E2[0:t0_jour] <- tendance_E2_1
data_journaliere$tendance_lin_E2[t0_jour:n] <- tendance_E2_2
ggplot(data_journaliere, aes(date, E2)) +
  geom_point(col = "grey") +
  geom_line(aes(y = tendance_lin_E2), col = "red") +
  geom_line(aes(x = date[(t0_jour)]), col = "blue") +
  ggtitle("Tendance linéaire avant et après l'ouverture de la 4e intersection")

Intersection 4

t0_jour <- which(data_journaliere$date=="2017-01-01")
N <- nrow(data_journaliere)
data_I4 <- data_journaliere[t0_jour:N,]
ggplot(data_I4,aes(x=date,y=E4)) +
  geom_point()

Détecter la cassure :

n<-nrow(data_I4)
chow_pvalues<-rep(NA,n)
for(t0 in 2:(n-2)){
 x1<-matrix(data_I4$date[1:t0],ncol=1)
 y1<-data_I4$E4[1:t0]
 
 x2<-matrix(data_I4$date[-(1:t0)],ncol=1)
 y2<-data_I4$E4[-(1:t0)]
 
 res_chow <-chow.test(y1,x1,y2,x2)
 chow_pvalues[t0]<-res_chow[4]
}

cassure <- which.min(log(chow_pvalues))
ggplot(data_I4, aes(date, E4)) +
  geom_point() +
  geom_line(aes(x=data_I4$date[cassure]), col = "blue")

Modèle linéaire sur les deux parties :

n <- nrow(data_I4)
data_I4_avant <- data_I4[0:cassure,]
data_I4_apres <- data_I4[(cassure+1):n,]

mettons en place les deux modèles

Valeur dans ‘data’ :

data_I4$tend_E4 <- NA
data_I4$tend_E4[0:cassure] <- data_I4_avant$tend_E4
data_I4$tend_E4[(cassure+1):n] <- data_I4_apres$tend_E4

Visualisation

ggplot(data_I4, aes(date, E4)) +
  geom_point(col = "grey") +
  geom_line(aes(y = data_I4$tend_E4), col = "red") +
  geom_line(aes(x=data_I4$date[cassure]), col = "blue") +
  ggtitle("Tendance linéaire avant et après l'ouverture de la 4e intersection")
## Warning: Use of `data_I4$tend_E4` is discouraged.
## ℹ Use `tend_E4` instead.
## Warning: Use of `data_I4$date` is discouraged.
## ℹ Use `date` instead.

On voit bien la différence ici. On observe une augmentation dans les premiers mois après l’ouverture de la 4e intersection, puis l’augmentation est beaucoup plus lente.

Ratios par intersection par jour et par mois :

Par jour, ratios entrées sur intersection

data_journaliere$E1_I1 <- data_journaliere$E1 / data_journaliere$I1
data_journaliere$E2_I2 <- data_journaliere$E2 / data_journaliere$I2
data_journaliere$E3_I3 <- data_journaliere$E3 / data_journaliere$I3
data_journaliere$E4_I4 <- data_journaliere$E4 / data_journaliere$I4

Par mois

data_mois$E1_I1=data_mois$E1/data_mois$I1
data_mois$E2_I2=data_mois$E2/data_mois$I2
data_mois$E3_I3=data_mois$E3/data_mois$I3
data_mois$E4_I4=data_mois$E4/data_mois$I4

Visualisation des ratios par mois :